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Abstract. Multiple scattering of polarised electromagnetic waves in diffusive media is in- 
vestigated by means of radiative transfer theory. This approach amounts to summing the 
ladder diagrams for the diffuse reflected or transmitted intensity, or the cyclical ones for the 
cone of enhanced backscattering. The method becomes exact in several situations of inter- 
est, such as a thick-slab experiment (slab thickness L 3> mean free path I 3> wavelength 
A). The present study is restricted to Rayleigh scattering. It incorporates in a natural way 
the dependence on the incident and detected polarisations, and takes full account of the 
internal reflections at the boundaries of the sample, due to the possible mismatch between 
the mean optical index n of the medium and that ri\ of the surroundings. This work 
does not rely on the diffusion approximation. It therefore correctly describes radiation in 
the skin layers, where a crossover takes place between free and diffusive propagation, and 
vice- versa. Quantities of interest, such as the polarisation-dependent, angle- resolved mean 
diffuse intensity in reflection and in transmission, and the shape of the cone of enhanced 
backscattering, are predicted in terms of solutions to Schwarzschild-Milne equations. The 
latter are obtained analytically, both in the absence of internal reflections (n = ni), and 
in the regime of a large index mismatch [njn\ C 1 or > 1). 
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1. INTRODUCTION 



Light undergoes multiple scattering when propagating through inhomogeneous media 
over distances much larger than one mean free path I. This may occur in a wide variety 
of situations, ranging from the atmospheres of stars and planets to biological tissues. The 
theory of multiple scattering of electromagnetic waves is an old classical area of physics 
[1-4], which has been developed for almost one century, mostly by astrophysicists. This 
subject has been experiencing an important revival of theoretical and experimental activity 
for one decade, motivated by the analogy between the effects of random disorder (weak 
or strong localisation) on the propagation of classical waves (electromagnetic, acoustic, 
seismic) and of quantum- mechanical waves (electrons in solids) . The first weak- localisation 
effect to be discovered has been the celebrated enhanced backscattering phenomenon, which 
takes place in a narrow angular cone around the direction of exact backscattering [5]. 

Typical laboratory experiments on multiple light scattering involve suspensions of 
polystyrene spheres or of Ti0 2 (white paint) grains in fluids. In these situations the mean 
free path I is usually much larger than the wavelength A of light, and the samples are 
often optically thick slabs, of thickness L ^> £. The regime of interest, i.e., A <C f < L, 
is characterised by a diffusive transport of radiation through multiple scattering. This 
diffusive regime admits three different levels of theoretical description, (i) The crudest 
approach is the diffusion approximation, where the multiply scattered intensity is described 
by means of an effective diffusion equation. The latter is only valid on length scales larger 
than £, so that is has to be supplemented by boundary conditions. As a consequence, this 
approach somehow keeps a phenomenological character, (ii) The mesoscopic approach, 
known as radiative transfer theory (RTT), has been used for long by the community of 
astrophysicists [1-4] . It is based on a local balance equation, keeping track of the direction 
of propagation of the intensity, (iii) The systematic microscopic approach consists in 
expanding the solution of the Maxwell equations in the random medium as a diagrammatic 
Born (multiple-scattering) series. 

In multiple light scattering experiments the quantities of most interest are the mean 
diffuse reflected and transmitted intensity, and the shape of the peak of enhanced backscat- 
tering. For these observables the diagrammatic approach greatly simplifies in the regime 
Acf <I. As far as mean quantities, averaged over the random positions of the scatter- 
ers, are concerned, the diffuse radiation is described by the sum of the ladder diagrams, 
which, in turn, amounts to RTT; the enhanced backscattering phenomenon is described by 
the so-called cyclical or maximally-crossed diagrams, which can also be summed up by an 
adaptation of RTT [6, 7]. On the other hand, the validity of vector RTT has been estab- 
lished on a rigorous basis [8], starting from a perturbative treatment of Maxwell's equa- 
tions, extending earlier developments [9] on multiple scattering of electromagnetic waves 
in plasmas. In the regime where the random fluctuations of the dielectric constant have 
short-range correlations, this approach rigorously justifies the use of the Schwarzschild- 
Milne equation of vector RTT with the Rayleigh phase function, which will be the purpose 
of the present work. 

The principles of vector RTT for electromagnetic waves, taking into account their 
polarisations, are exposed in the book by Chandrasekhar [1] , which also contains a formal 
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analytical derivation of the diffuse intensity for Rayleigh scattering, in the absence of 
internal reflections. This approach is needed in order to obtain predictions at a quantitative 
level, concerning observables like the diffuse intensity in reflection and in transmission, and 
the enhanced backscattering cone. In particular the diffusion approximation alone cannot 
yield such accurate predictions, chiefly because boundary conditions cannot be dealt with in 
a fully satisfactory way. Surprisingly enough, in the modern era of multiple scattering only 
very few authors have used RTT. Thus far the major investigations of the weak localisation 
of light, including polarisation effects, have rather used either the diffusion approximation 
[10-12] or numerical simulations [13]. Several other bulk properties of multiple scattering 
of electromagnetic waves have been investigated along these lines, including especially the 
effects of Faraday rotation [12, 14, 15] and of absorption [16]. Exact results on polarisation 
effects on the backscattering cone have only appeared very recently. Mishchenko [17] has 
derived general properties of the behaviour of polarisations under time reversal, obtaining 
thus for the first time a consistent derivation of the enhancement factors in the direction of 
exact backscattering. The full shape of the backscattering cone has then been investigated 
by Ozrin [18], who did not, however, come up with a full analytical solution of the latter 
problem. In previous works, we have considered the case of scalar waves undergoing 
multiple isotropic [19, 20] and arbitrary anisotropic scattering [21]. We have shown how 
RTT takes proper account of the skin layers, where light is converted over a few mean 
free paths from a free beam to a diffusive field and vice-versa, and how it allows to deal 
with the effects of internal reflections due to the index mismatch at the boundaries of the 
sample. The latter effect has been the subject of much activity recently [22-26] . 

The goal of the present paper is to extend our investigations to the multiple scattering 
of electromagnetic waves, obtaining thus for the first time a complete analytic description 
of the diffuse intensity and of the backscattering cone in the regime A C ^ < I, including 
both polarisation effects and internal reflections. We shall restrict the analysis to Rayleigh 
scatterers for definiteness. Section 2 contains general results on vector RTT. The observ- 
ables of interest, with their dependence on polarisations and index mismatch, are expressed 
in terms of solutions of appropriate Schwarzschild-Milne equations. Sum rules and other 
general properties are given. These predictions are then made more quantitative in two 
situations: (a) in the absence of internal reflections (section 3), where a full analytical 
solution of the problem is given; the results of refs. [1, 17, 18] are made more precise and 
put in a broader perspective; (b) in the opposite regime of a large mismatch of optical 
index between the sample and the surroundings (section 4). Section 5 contains a brief 
discussion of our findings. 

2. GENERALITIES 

2.1. General formalism 

Throughout the following we consider the multiple scattering of electromagnetic waves 
by a diffusive medium containing a low density p of identical Rayleigh scatterers charac- 
terised by their cross-section cr, so that the mean free path I = l/(pa) is much larger than 
the wavelength A of radiation in the medium. The diffusive medium has the form of a slab 
(0 < z < L), infinite in the transverse directions. We introduce the optical thickness b of 
the sample through L = b£, and the optical depth r of a point in the sample (0 < r < b) 
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through z = t£. We shall consider either optically thick slabs (6 > 1), or semi-infinite 
samples (b = +00). We investigate the general situation where the mean optical index n 
of the sample is different from that n\ of the surrounding medium. Whenever there is an 
index mismatch (m = njn\ 7^ 1), internal reflections take place at the boundaries of the 
sample. Useful definitions and notations are summarised in Table 1. 

We closely follow the definitions and notations of Chandrasekhar [1] . We measure the 
polarisation of the radiation in a fixed frame, introducing spherical co-ordinates (see Table 
1). For radiation propagating in the angular direction (6>, tp) with respect to the z-axis, 
the complex components of the electric field E in the plane transversal to that direction 
will be denoted by (Eg.E^). The component Eg is parallel to the unit vector and 
contained in the meridian plane, defined by the direction of propagation and the normal 
to the boundaries of the sample, while the component E v is parallel to the unit vector <p> 
and normal to the meridian plane. The alternative notations (E\\,E±) and (Ee,E r ) can 
be found in the literature. We introduce the following vector of four Stokes parameters, 
or Stokes vector for short 

11 = \Eg\ 

12 = l^vl 2 / 91 v 

1 3 = U = 2Re (EgE^) ^ L> 

1 4 = V = 2lm (E E;). 

Here and throughout the following, boldface symbols represent 4-component vectors or 
4x4 matrices. The description of polarised radiation by means of Stokes parameters is 
very commonly used [27]. This formalism has many advantages: the Stokes parameters 
add up for light beams superposed incoherently; a scattering event is described by a linear 
transformation of the Stokes parameters, i.e., by the action of a scattering matrix on the 
Stokes vector I. We finally recall the definition [27] of the degree of polarisation P of 
radiation described by a Stokes vector I: 



II + I2 



2.2. Schwarzschild-Milne equation 

For reasons exposed in the Introduction, we shall use RTT to investigate the average 
reflected or transmitted intensity in the regime £ 3> A. We consider first the situation of a 
semi-infinite medium, for simplicity. 

The mean diffuse radiation propagating in the direction (6, ip) in the medium at depth 
t = z/£ is described by its Stokes vector I(r, fx, <p), with the notations (see Table 1) 

/ u = cos6>, v = sin6> = a/1 — \j? . (2.3) 
Along the lines of ref. [1] , the RTT equation reads 

d 



4 



where the vector source function T(r, fi, ip) is defined as 

r(r, /*, = J 1 ^ jf* ^P(^ J, <p').I(T, *,'), (2.5) 

with the Rayleigh phase matrix P(fi, <p, fi', <p') being the matrix describing a scattering 
event, expressed in the fixed frame related to the sample. Its explicit form will be given in 
eqs. (2.10), (2.11). 

The RTT equation (2.4), with appropriate boundary conditions, leads to the following 
linear integral equation for the source function, referred to as the Schwarzschild-Milne (SM) 
equation 

r(r, ^, (p) = P(a», V, A*a, ^ a ).Iae" T/Ma 



/ /-27T J / 



The right-hand side of this equation has the following interpretation: 



(2.6) 



• The first line is the exponentially damped contribution, with a suitable normalisation 
[19, 20], of the incident beam, characterised by an incident direction (9 a ,<p a ) and a Stokes 
vector I a ; 

• The second (third) line is the bulk contribution of diffuse light scattered from a smaller 
(larger) depth r'\ 

• The fourth line is the layer contribution of diffuse light scattered from depth r' and then 
reflected at the boundary (r = 0). The effect of the boundary is described by a reflection 
matrix R(/(i) and a transmission matrix T(/u), namely 



R(/0 



( hi(^)l 




V o 









V±(p)\ 2 o 

Re(r||(At)r ± (//)*) - Im(r|| (//)r ± (//)*) 

Im(r||(At)r ± (//)*) Re (ry (fj)r ± (//)*) / 



T(/*) = 



m 2 z/ 2 






V o 





\t±<J*)t 
















Re(t||(/i)t ± (/i)*) 
Im(t||(//)t ± (//)*) 



-Im(t||(/x)t ± (/i)*) 
Re(*ll(//)t ± (//)*) 



/ 



(2.7) 
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In these expressions r^(n), r±(fjL) and t±(n) are the Fresnel reflection and transmis- 

sion amplitude coefficients, respectively. The latter only depend on the inner incidence 
angle 9 and on the index mismatch m, according to 



my/1 — m 2 i> 2 



\i + my/1 — m 2 v 2 ' 

2V1 - m 2 v 2 
\i + my/1 — m 2 v 2 ' 



r±(lJ>) = 



VT 



m 2 v 2 



mLi 



y/1 — m 2 v 2 + m/i 
2y/l-m 2 v 2 



(2.8) 



y/1 — m 2 v 2 + rail 



In the case of partial reflection (see Table 1), these coefficients are real, with absolute 
values less than unity. In the case of total reflection, the reflection coefficients are pure 
phases, i.e., complex numbers with unit modulus, while the transmission coefficients vanish 
by convention. The first two diagonal elements of the reflection matrix R(/i) of eq. (2.7) 
read 

|r||(^)| 2 = R^li) = 1 -T||(/u), |r ± (^)| 2 = R ± (ii) = 1 - T ± 0*), (2.9) 

in terms of the Fresnel reflection and transmission intensity coefficients. 

It is advantageous to expand the 99-dependence of all quantities in the complex trigono- 
metric polynomials {e lfc(,p }, with — 2 < k < 2. We thus set 



I(/W) = lW We^, Tfatp) = r (fc) (^)e^, 

fc=-2 fc= 
2 



fc=-2 



fe=-2 



The Rayleigh phase matrices read 



P ( °W) = 7 



/2(l-/x 2 )(l-// 2 ) + /ry 2 /« 2 \ 

Li' 2 10 0' 



V 2/jl/jl'J 



p( 1 )(^, / i') = P ( - 1) *(/*^ / ) = V 



/ z/i 0\ 



-2i/x' 1 

V 001/ 



/ ii 2 ii' 2 -fi 2 ifj 2 ^ 0\ 



P< 2 W) = P<- 2 >W) = i 



,'2 



fj,'" 1 —ill' 
2z / u / u' 2 2i/j 2fjL/j,' 

V 000/ 



(2.10) 



(2.11) 



and the SM equation (2.6) splits into the following five decoupled integral equations (—2 < 
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k < 2) 

r <*> (r, n) = P« fa, fi a ).I a e- tk ^ ~ T/ ^ 

+ f dr ' f 1 ^" (TV)/ ^ w (M').r w (r>') 

Jo Jo ^A* 



+ / °° dr/ / ^7^ (T '- r)/M 'P (fc) (^V).r( fe )(r',V) 

+ / 1 ^-^')^R( / .').pW(^').r W (r^'). 
Jo Jo -^A* 



2.3. Solutions to the SM equation and sum rules 

We now turn to the analysis of the solutions to the SM equation (2.6). We shall 
investigate their general symmetry properties, and show that they obey some remarkable 
sum rules. 

We start by introducing the matrix Green's function G(r, /i, <p, r', fi', <p'), defined as 
being the solution, which remains bounded as either r or r' goes to infinity, of the SM 
equation with a matrix ^-function source term, namely 

G(t, fx, if, t\ if') = P(jm, if, fjf, (p')S(t - t') 
+ I dr" 



Jo dT "l V 7 / ^e-^-^^^^^^O-GC^^^ 
+ 1^ dr" jf 1 ^ fj ^ e -(-"-)/M"p( Mj v,, V', ¥>").G(r", V', ¥>", r', <p') 

+ ir dT " Jo W r ^ e " (r+r " )/M " R (^)- P (^' V. A<", v")-G(r", p", r', p', </). 

(2.13) 

The source function T(r, /i, </?), solution of the SM equation (2.6), is then given by 

V, <P) = I*, V 7 , l*a, (p a ).Ia, (2-14) 

where 

I-+00 i 

r(T,fJL,(p,fJL a ,<p a ) = dr / e" r /Ma G(r, n,(p,r', n a ,<Pa)- (2.15) 

Jo 

We also define the matrix of bistatic coefficients, or bistatic matrix for short 



r+oo 

ry(Ha,<Pa,Vb, <Pb) = / dre~ r/At! T(T, -fJL b ,(p bt fl at (p a ) 
JO 



r+oo r+oo 

/ dre"^* / dr'e^'/^G^,-^,^,/,^,^). 
Jo Jo 



(2.16) 



The invariance of the Rayleigh scattering mechanism under time reversal implies the 
following symmetry properties of the quantities defined so far. We introduce the constant 
matrices 

/l 0\ /l 0\ 



K = 







1 


Vo o 6 



o o \ 



L = 



/ 



1 


Vo o 




-i 



(2.17) 







the matrix K being denoted by Q 1 in ref. [1]. The Rayleigh phase matrix P(p, (p, n' , ip') 
has the symmetry property = 1, • • • , 4) 



(K.P)ij(ji, (p, fjf, if') = (K.P)^^', ip', /i, ip), 
(L.P)ij(ji, (p, //, ip') = (L.P)^(-//, if', -fi, if). 



(2.18) 



It then follows from their definitions (2.13), (2.16) that the matrix Green's function and 
the bistatic matrix obey the symmetry relations (i, j = 1, • • • , 4) 



(K.G)y(r, ii, if, r', n' , if') = (K.G)ji(r', <p', r, <p), 
(L.G)y (r, f/, </?, r', p') = (L.G)^(r', -//, <//, r, -/x, p), 
(L.7) y (ji, ip, v?') = (L. 7 )^(//, ip', n, <p). 



(2.19) 



We now investigate the asymptotic behaviour of quantities deep inside the medium, 
i.e., for r — > +oo. It is expected on physical grounds that the diffusive medium depolarises 
the incident radiation, so that both \{r,n,(p) and T(r,n,ip) become proportional to the 
Stokes vector of natural (unpolarised) light, namely [27] 



Inat — 



1 



Vo/ 



(2.20) 



This assertion will be made quantitative in section 3, where the extinction lengths of the 
all other modes will be determined. The asymptotic behaviour of the matrix solution 
r(r, ii, ip, fx a , tp a) then assumes the form 



T{t,h,<p, li a ,<fa) 



/n(A»o) 72 (aO o o^ 

Tl(AO ^(aO 



V o ooo/ 



(r — > +oo). 



(2.21) 



Furthermore it should be noticed that the homogeneous SM equation (2.6), (2.12), without 
a source term, has a vector solution Th(t, h) in the (^-independent (k = 0) sector, growing 
linearly as r -> +oo. We shall refer to the latter solution, normalised as 



r#(T, fj) « (t + Tq) I nat , 



(2.22) 
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as the homogeneous solution, for short. 

The constant r and the functions ti(ji) and r 2 (//), which show up in eqs. (2.22) and 
(2.21), respectively, are unknown so far. It will be shown that they bear the full non- 
trivial dependence of the physical observables on the index mismatch. They will also be 
determined analytically in the absence of internal reflections (section 3) and in the large 
index mismatch regime (section 4). 

For the time being, we pursue our investigation of general properties. The special 
and homogeneous solutions to the SM equation in the k = sector can be related among 
themselves as follows. The column vectors of the matrix Green's function obeying eq. 
(2.13) become proportional to the homogeneous solution, as either r or r' goes to infinity: 

^im G^(t,^t'^') = ^-(T h Ut^) (z,j = 1,---,4), (2.23) 

where the constant D will be determined and interpreted in a while. As a consequence of 
the above definitions, we have 

tM= lim U > l ^ =— / dre-^(rH)i(r, -/*») (<,j = l,2). (2.24) 

We end up by deriving two groups of sum rules obeyed by the quantities defined above, 
which are related to the F- and X-integrals, with the notations of ref. [1]. 

• First, we consider the F-integral, defined as 

F(r) = ^(l^^^+Ur^^)). (2.25) 

The RTT equation (2.4) implies dF/dr = 0, expressing thus the conservation of the flux 
in the ^-direction. We consider first the F-integrals F\ and F2 associated with the special 
solutions Tii and T^, i.e., the first two column vectors of the matrix (2.15). The r — > +00 
limit determines F\ = Fi = 0; the r = values then yield the sum rules 



1 d/i 
d^ 



T \\ (^)7ll ^a)+T 1 _ {jl,Ha) 



(2.26) 



The n a — > +00 limit of these equations, using eq. (2.24), yields 

11 d/i 



J ^[mn)rM + T ± ^)rM]=l. (2.27) 



Similarly, we consider the F-integral F H associated with the homogeneous solution T H . 
This does not yield any independent sum rule, but rather leads to the determination of 
the unknown constant D, namely 

D = \, (2.28) 
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to be interpreted as the dimensionless diffusion constant, i.e., -D p h ys = c£/3 in physical 
units. 



• Second, we consider the FJ-integral, defined as 

- 1 d^ ^ 



K(r) = J i J J g (l x (r, fx, <p) + I 2 (r, /*, p)) . 



(2.29) 



The RTT equation (2.4) implies dif/dr = — F, whence if(r) = — Fr + F"o, with if being 
a constant. Along the lines of the above derivation, we obtain the sum rules 



L 



1 d/x 

T' 



(1 + i?||Gu)) 7 ^V /O + (1 + iU/^hKV //«)J = gTiC/Xa) + ^, 



(2.30) 



(l + -R||(/*))721 V>/V) + (1 + -R±(^))722 (^»^a)J = gT 2 (^a) + 



and 



1 d/x 



(l + J R||(/x))r 1 (/x) + (l + J R ± (/i))r 2 (/x) 



= T . 



(2.31) 



The sum rules (2.30), (2.31) express that the multiple-scattering problem in a semi-infinite 
sample is invariant if a finite slab of any thickness is added, or removed, from the sample 

[!]• 

2.4. Diffuse reflected intensity 

The diffuse reflected intensity for a semi-infinite sample can now be calculated, along 
the lines of refs. [19-21]. The incident radiation is characterised by the direction (9 a ,^ a ) 
and the Stokes vector I a ; the reflected radiation is detected in the direction (9b, <fib) and 
in a polarisation state characterised by the Stokes vector I5. Our prediction for the mean 
reflected intensity per solid-angle element reads 

dR(a — > b) A R, n n cos 9 a , , . , , 

— = A (9 a ,ip a ,9 b ,ipb) = 47rm 2 (Ib|T(/Xfe).L.7(/i a , (f a , f^b, <Pb)-T(fta) |Ia)- 

(2.32) 

In the absence of index mismatch we obtain the simpler expression 

A R (9 a ,(p a ,9 b ,(pb) = -r~ (l6|L.7(/Lt a ,^ a ,// b ,y? 6 )|l a ). (2.33) 



2.5. Diffuse transmitted intensity 

The diffuse transmitted intensity through an optically thick slab (& > 1) can also be 
calculated along the lines of refs. [19-21]. A first step consists in building up the solution 
T(b, r, fi, <p, fi a , <fa) of the SM equation pertaining to the thick-slab geometry. This solution 
can be expressed in terms of the solutions T(r, fi, <p, fi a , cp a ) and Th(t, fx) pertaining to the 
semi-infinite geometry, by means of a matching procedure. It turns out that only the (1,2) 
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sector of the matrix solution matters, since all the other matrix elements are exponentially 
small in the optical thickness b. We thus get = 1, 2) 



Tij(b, T,H,ip,fJ, a , (fa) 



T, 



b + 2r 



(r H )i(T,//) (r finite, b-r > 1), 

(b-r finite, r > 1). 

(2.34) 

Both expressions lead to a linear (diffusive) behaviour in the bulk of the sample (r ^> 1, 
b — t 3> 1), namely 



T«(//q) 

6 + 2r 



1^ (6, r,//, /ia^a! 



6 + r - r 
6 + 2r 



(i,j = l,2). 



(2.35) 



The mean diffuse transmitted intensity through an optically thick slab can then be 
derived explicitly, again along the lines of refs. [19-21]. The incident radiation is charac- 
terised by the direction (0 a ,(p a ) and the Stokes vector I a . The transmitted radiation is 
detected in the direction (0 b , ip b ) and in a polarisation state characterised by the Stokes 
vector lb. Our prediction for the mean transmitted intensity per solid-angle element reads 



dT(a^b) A T (0 a ,0 b ) 



dttt 



b + 2r ' 



(2.36) 



with 



A T (0 a ,0 b ) = 



COS 6 n 



(l 6 |A(/i a ,// 6 )|l a >. 



12nm 2 Hafih 

The matrix A(// a , fj, b ) has non-zero elements only in the (1, 2)-sector, namely 



(2.37) 



A(ji a ,ji b ) 



( T^{P'a)ri{p,a)T\\{ii b )T 1 {ji b ) T||(^ a )ri(^ a )T ± (^ 6 )r 2 ( / u 6 ) 
rj_(/Ua)r 2 (/U a )T||(^ b )ri(^ b ) Tj_(^ a )r 2 ( / u a )T ± (^ b )r 2 ( / u 6 ) 



V 



o 
o 
















0/ 



(2.38) 



2.6. Enhanced backscattering cone 
2.6.1. Generalities 

In the regime A £ of interest, the enhanced backscattering phenomenon takes place 
in a narrow cone around the exact backscattering direction, of angular width of order X/£. 
As recalled in the Introduction, the shape of the cone of enhanced backscattering for a 
semi-infinite medium is given by the sum of the cyclical, or maximally-crossed, diagrams. 
This summation can be performed by means of an adaptation of RTT. This property has 
been exploited extensively in the case of scalar waves [6, 7, 19-21]; it has been extended 
more recently to polarisation effects for electromagnetic waves [17, 18]. 
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We restrict ourselves to a semi-infinite medium and to normal incidence (9 a = 0). We 
define the dimensionless transverse wavevector of the outgoing radiation as 

Q = q£, (2.39) 

with a magnitude 

Q = q£ = k£6 = k x l9 u (2.40) 

with 9i being the observation angle. We assume for definiteness that the vector Q is par- 
allel to the x-axis, namely Q = Q% with Q > 0. In order to cure the ill-definedness of 
the co-ordinate system at strictly normal incidence, we choose to give the initial wavevec- 
tor an infinitesimally small positive component along x. We thus set 9 a = + , (p a = 0, 
so that 6 = x and up = y. We then introduce a Q-dependent matrix of bistatic coeffi- 
cients, jij(Q, Ha, <fa, A*6, (£>b). The latter is defined, in analogy with eq. (2.16), in terms of 
the matrix source function I\j (Q, \x, <p, [i ai <p a ). This matrix solves the Q-dependent SM 
equation, obtained by replacing in eq. (2.13) the exponential damping factor exp(— r/fi') 
by exp ( — (1 — iQ.n)r/ / u / ), where n is the unit vector in the direction (6>',<//), so that 
Q.n = Qv' cos ip' . 

We now turn to the explicit shape of the enhanced backscattering cone. It can be 
expressed [17, 18] in terms of the values at normal incidence of the bistatic coefficients, 
lijiQ) = lijiQi^a = l>V?a = 0, [ib = l,<fib = 0). To be more specific, the total reflected 
intensity near the backscattering direction, i.e., for 9 <^ 1, kl ^> 1, and Q = k£9 > fixed, 
reads 

A(Q) =A L + A C (Q) - A ss = — _L_(I 6 |L.( 7 L + 1 C (Q) - l SS )\l a ). (2.41) 



• The first term in eq. (2.41), given by the sum of the ladder diagrams, coincides with the 
expression (2.32) for the background reflected intensity. At normal incidence it assumes 
the general form 



1 L = li^a = 1, <Pa = 0, (lb = 1, Vb = 0) = 



/7n 


712 





° \ 


712 


7n 














712 — 7n 





V 








744/ 



(2.42) 



where the three constants 711, 712, and 744 only depend on the index mismatch. 

• The second term in eq. (2.41) is given by the sum of the maximally crossed, or cyclical, 
diagrams. It represents the contributions of the interference between the sequences of any 
number (N > 1) of scattering events and their time-reversed counterparts. At normal 
incidence we have [17, 18] 



7 C (Q) = 



/7n(<2) 712(G) 

712 (Q) 722 (Q) 


V 






733 (Q) 









744(G) J 



(2.43) 
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with 



7i2(<2)= g (744(g) -73 3 (Q)), 

733(Q) = \ (733(Q) +744(Q)) -7i2(Q), 

744(g) = ^ (733(g) + 744(g)) + 712(g). 



(2.44) 



• The subtracted third term in eq. (2.41) is the contribution of the single-scattering 
(N = 1) events, which are their own time- reversed counterparts, and must not be double- 
counted. At normal incidence it reads 



7 SS = 



/l \ 

10 

0-10 

Vo -1/ 



(2.45) 



The actual calculation of the g-dependent bistatic matrix j(Q, fj, a , <p a , Hb, ^Pb) goes as 
follows. By expanding the g-dependent SM equation in the trigonometric polynomials 
{e lk(p }, we obtain the following system of coupled equations (—2 < k < 2) 

r«(r,^) = P«(^,^ a ).I a e-^-^ 

„ r „i 2 



Jo Jo Z A* J= _ 2 
+ Hdr' r^e-^'-^'P^^, -//)• E i*- J "J*- i (Q(r'-r)i// / *')r«)(r / , V) 

Jr JO Z A* j = -2 



+ /"dr'/^e-^^'R^.P^^^)- E i*-^*-i(Q(r + r / )i///i / )r«)(r / ,/i'), 
Jo Jo „~; 9 



j = -2 



(2.46) 

where z/ has been defined in eq. (2.3), and where the J n (z) are the Bessel functions, which 
admit the integral representation 



i J n (z) = / — — exp (zz cos ip — tn(pj , 

JO 27T 



(2.47) 



and possess the symmetry property 

J_ n (z) = J n (-z) = (-l) n J n (z). (2.48) 

2.6.2. Linear polarisations 

We now investigate the case where the initial beam is linearly polarised, and a linear 
polarisation of the outgoing beam is detected. Let ip a and ipb be the respective angles 
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between the directions of the polarisations and the direction of the Q- vector, i.e., the 
positive x-axis. The corresponding Stokes vectors read 



f COS 2 V>a \ 

sin 2 ip a 
sin(2^ a ) 

V o / 



I COS 2 -06 \ 

sin 2 ip b 
sin(2ip b ) 

V o / 



(2.49) 



By inserting these expressions into the results (2.41-45), we obtain that A L and A ss only 
depend on the relative angle 



between the directions of both polarisations, namely 

4 



A L 

A ss 



Tr(m + 1) 
3 



(711 cos 2 * + 712 sin 2 *) , 
cos 2 



n(m + l) 4 

whereas A C (Q) depends separately on both polarisation directions: 
4 



(2.50) 



(2.51) 



7n (Q) cos ip a cos V6 + 722(Q)sin VaSin 1p b 
+ (2^f 12 (Q) - ^3s(Q) - 744(Q)) cos ifj a sin ifj a cos ifj b sin ifj b (2.52) 
+ o (744 (<3) - 733 (Q)) (sin 2 ip a cos 2 ip b + sin 2 ip h cos 2 ip a ) 



We define as usual the enhancement factor B(Q) as the ratio between the total re- 
flected intensity and its background value: 



B(Q) 



A L + A C (Q) - A ss 
A 1 ' 



(2.53) 



Right at the top of the backscattering cone, corresponding to the exact backscattering 
direction {Q = 0), the expressions (2.41), (2.43) simplify to 



A C (0) = \ ( 7ll - \ ( 7ll + 7l2 - 744) sin 2 fr) . 
7r (m + 1 r \ 2 J 



:(m+ 1) 

The enhancement factor thus reads 



5(0) 



2711 + \ (- 3 7ii + 7i2 + 744 + ^ ) sin 2 # 

7n + (712 - 7n) sin 2 ^ 
14 



(2.54) 



(2.55) 



The maximal enhancement factor B\\ is observed for parallel detection, i.e., ^> = 0, whereas 
the minimum B± corresponds to perpendicular detection, i.e., \1/ = ±7r/2. These extremal 
values read 



Bn =2 



4 7 



5± = 



11 



711 + 712 + 744 

2712 



(2.56) 



A celebrated and universal feature of the enhanced backscattering cone is the tri- 
angular shape of its top. Within the present formalism, and in analogy with previous 
studies [19-21], this phenomenon is described as follows. For Q< 1, and for i,j = 1,2, 

the solution T^(Q, r, fx) of the Q-dependent SM equation has a term linear in Q that is 



proportional to the homogeneous solution (Th)i(t, /u), namely 

T^Q, r, fj) = Ttj (r, fi) - C % Q(r H )i(r, fi) + 0(Q 2 ) = 1, 2). 



(2.57) 



The constants C; are then fixed by requiring that the above solution falls off as exp(— Qt) 
for Qt 3> 1. This general property will be checked explicitly in section 3 in the absence of 
internal reflections. We thus obtain C\ = C 2 = Ti(l) = r 2 (l), so that 



Hj(Q) = Hj ~ -r^l) 2 Q + 0(Q 2 ) = 1,2), 



(2.58) 



and finally 

A C (Q) = 



(711 " \ (711 + 722 - 744) sin 2 * - \ti(X) 2 cos 2 * Q + 0(Q 2 )^j . 

(2.59) 



7r(m + 

Along the lines of refs. [19-21], we define the width AQ of the triangular cone as 



Q 



A^(Q) = A^(0)[l--^ + om 



(2.60) 



The sharpest cone, namely the smallest width AQ, is observed for parallel detection, i.e., 
\1/ = 0, where we have 

AQ, = J^p. (2.61) 

The universal features of the top of the enhanced backscattering cone described so far 
only depend on Q and The full shape of the enhancement factor B(Q) weakly depends 
separately on the directions ifj a and tpb of both polarisations. This phenomenon will be 
illustrated in section 3.2 in the absence of internal reflections. 

2.6.3. Circular polarisations 

We end up by investigating the case of circularly polarised beams at normal incidence. 
The corresponding Stokes vectors now read 



(2.62) 









2 



, I& = 


2 



w 




W 
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where the helicity is a a = 1 (respectively, a a = — 1) if the incident beam has a left 
(respectively, right) circular polarisation, and similarly for the helicity of the detection 
channel. By inserting these expressions into the results (2.41-45), we observe that the 
various backscattered amplitudes only depend on the relative helicity 



according to 



a a a b , (2.63) 



A ss = , - 

7r(m + l) 4V ; ' 

7r(m + l) 4 (264) 
= 7 1 ru hn(Q) + 722(Q) - 733(Q) + 74 4 (Q) 



7r(m + l) 4 



+ (2 7 i2(Q) +73 3 (Q)+74 4 (Q))S 



The enhancement factor B^(Q), defined in analogy with eq. (2.53), is larger for the 
helicity-preserving channel (E = 1) than in the channel of opposite helicity (E = — 1). In 
particular, right at the top of the cone, we have 

B 1 = 2, j?_ 1 = 37ll ~ 7l2 ~ 744 ~ f . (2.65) 
711 + 7i2 - 744 

The maximal value of the enhancement factor in the helicity-preserving channel is exactly 
equal to two, because A L = A c (0) and the single-scattering contribution vanishes. Cor- 
rections to this exact factor of two for denser diffusive media (£/X not very large) have 
been measured in a recent experiment [28] , and given a theoretical interpretation in terms 
of recurrent double scattering [29] . 

Another consequence of the result (2.64) is that the characteristic triangular shape of 
the cone only shows up in the helicity-preserving channel. The associated width, defined 
in analogy with eq. (2.60), reads 

AD 3(711 +712 + 744) /O R«^ 

= T^rp • (2-66) 



3. EXACT SOLUTION IN THE ABSENCE OF INTERNAL REFLECTIONS 

This section is devoted to the exact solution of the various SM equations introduced in 
section 2, in the case where there is no optical index mismatch between the sample and the 
surroundings, so that there are no internal reflections: the reflection matrix R(/u) vanishes. 
Therefore the SM equations (2.6), (2.12), (2.46) involve convolution kernels, which only 
depend on the difference of optical depths t — t'. The problem is, however, still non-trivial 
because of the semi-infinite geometry (0 < r < +00). We have found it worthwhile to 
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expose a self-contained derivation of the Wiener-Hopf technique, and of the results known 
previously, and already exposed in the book by Chandrasekhar [1]. 

The vector RTT problem is considered in section 3.1. The outcomes concerning diffuse 
reflection and transmission are compared in detail with those corresponding to multiple 
isotropic scattering of scalar waves [19, 20]. Section 3.2 deals with the enhanced backscat- 
tering phenomenon. We derive closed-form expressions for the five functions describing 
the full shape of the enhanced backscattering cone, up to the numerical solution of the 
9x9 system (3.74). The present analysis thus goes one step further than the recent work 
by Ozrin [18]. 

3.1. Diffuse reflection and transmission 

In this section we derive the exact solution to the SM equations (2.6), (2.12) in the 
absence of internal reflections, obtaining thus predictions for the diffuse reflected and 
transmitted intensity. We introduce the following parametrisation 



r(°W) 



/A(t) + B(t)(1-^)\ 

Mr) 


V C(r)/x J 

\ E(r) 



(3.1) 



( * \ 

-l 

—2ifx 

V o / 



where v has been defined in eq. (2.3). The functions A(r) 
equations 

3 



F(t) obey the integral 



A 
B 
C 
D 
E 
F 



((/^Ii + I 2 ) e~ T '^ + (M + M 2 ) * A + (M 2 - M 4 ) * b) , 
| (((2 - 3/12)1! - I 2 ) e - r /^ + (M - 3M 2 ) *A + (2M - 5M 2 + 3M 4 ) ★ s) , 

(^ a I 4 e- r/ ^+M 2 *c), 
| (z/ a (2^ a l! + zl 3 )e- r /^ + (M + M 2 - 2M 4 ) * £>) , 
= - A (z/ a I 4 e- T /^ + (M - M 2 ) * E) , 

= | ((^i! - I 2 + z^ a I 3 )e- r/ ^ + (M + 2M 2 + M 4 ) * f) , 



(3.2) 

where the brace shows that the equations for A(t) and B(r) are coupled, while the other 
four are decoupled. In the above equations, the star denotes the convolution between a 
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kernel M(t — r') and a function A(t), defined as 

r+oo 

{M*A){t)= M(t-t')A(t')<1t'. 
Jo 

The kernels entering eq. (3.2) are the following even functions 

M 2p (r) = f ^e-W". 
Jo 



(3.3) 



(3.4) 



3.1.1. Preliminaries 

As recalled above, the integral equations (3.2) are exactly solvable because of their 
convolution structure, which suggests to utilise the Laplace transformation. Along the 
lines of refs. [19-21], the Laplace transform of a function A(t) defined for < r < +oo 
will be denoted by a(s) (the corresponding lower-case letter), and defined as 



/■-t-oo 

a(s) = / A(r)e ST dr, 
Jo 

while the Laplace transform of the kernels M 2p (r) read 

m 2p (s) = / M 2p (r)e ST dr = dfi- ' 

J-oo Jo 1 



(3.5) 



9 9 



(3.6) 



i.e., explicitly, 



m (s) 



1 . 1 + s 
— In . 

2s 1 - s 



m 2 (s) = —(m (s) - l), 
s 

If s 2 
m 4(s) = im (s) - 1 - — 



(3.7) 



We also define for further use the following linear combinations of the above kernels 



01 (s 

<p2 (S 

4 (S 
05 (S 



1 - ^(mo(s) -m 2 (s)), 



1 



1 - g ( m o( s ) ~ m 2(«)) 



1 - -(m (s) + m 2 (s) - 2m 4 (s)), 
3 

1 - o ( m o( s ) + 2m 2 (s) + m 4 (s)), 

o 

1 - ^m 2 (s), 



(3.8) 
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which we shall refer to as the kernel functions. Both the kernels m 2p (s) and the kernel 
functions (j) n (s) are even functions of s, analytic in the s-plane cut along the real axis from 
— oo to —1 and from +1 to +00. 

In the following we shall need to factor the (j) n (s) into the corresponding so-called 
Wiener-Hopf .//-functions, defined after ref. [1] by the identity 

MS) = H n (s)H n (-s) (* = 1 '-' 5 )' ^ 

together with the condition that H n (s) is analytic in the left half-plane Re s < 0. Consider 
first the case of a rational function of the form 

M 



a=l 



N 



(3.10) 



IK* 2 -pD 



b=i 



with 2M zeros and 2N poles at arbitrary positions, with Rez a > 0, Kepb > 0. The 
factorisation (3.9) is elementary in this case, and the associated .ff -function reads 

N 

IK*-^) 

H(s) = ^ . (3.11) 



o=l 



This expression can be recast as a complex contour integral, yielding thus an explicit 
representation of the if-functions in the general case: 



ff n (a)=exp(/ ^ 



dz <f>' n (z) 



4>n{z) 



\n(z — s) ) = exp I — / — : 



dz In 4> n (z) 



2%i z — s 



(Res<0), (3.12) 



where the vertical contour can be placed at Re z = 0. In the present case it is advantageous, 
especially for the purpose of numerical evaluation, to change variables from z to an angle 
(3 such that z = i tan (3. We thus get 



with 





H n (s) = 


e*J S r /2 d8 ) 
\k Jo sin 2 13 + s 2 cos 2 13 ) 


(Re 


MP) 


4^ 


/ (cot 2 /3 + l)(/3cot/3-l) + l), 




MP) 


= cot 2 (3 


1 - ^ ((cot 2 /3 + 1)03 cot /3- 1) + l) 


> 


MP) 




'{2cot A (3 + cot 2 (3- \){(3cot(3- 1) + 


2 

-cot' 
3 


MP) 


—§( 


\cot 2 [3- l) 2 (/3cot/3- 1) + ^cot 2 /3 + l), 


MP) 


= 1 + - cot 2 /3(/3 cot /3 - 1). 





(3.13) 



P-l), 



(3.14) 
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The following values of the .ff-functions will play a role hereafter. First, the kernel 
functions have the following series expansions around the origin 



xi ^ 1 1 2 



,.1 3 o i / \ 3 13 o 

>(s) = — I s H , 0s(s) = s + 

5 35 > ^ ^ 10 70 



/ / \ " 2 / / \ 1 ^2 

^ 4(S)= 10-70 S + '"' ^ = 2" 10* + 
Eq. (3.9) yields H n (0) = l/y/(f> n (0), hence 



(3.15) 



#l(0) = A # 2 (0) = A ^3(0) = ff4(0) = \/y, ^ 5 (0) = V^. 




(3.16) 



Second, for large s, namely |s| — > +oo with Res < 0, the functions H n (s) with n ^ 2 go 
to unity, while we have H 2 (s) ~ — s. Finally, the values of the if-functions at s = — 1 can 
be accurately determined from the integral representation (3.13), (3.14). We thus obtain 



#i(-l) = 1.277973, # 2 (-l) = 3.469485, H 3 (-l) = 1.465 877, 
H 4 (-l) = 1.396 266, H 5 (-l) = 1.203 622. 



(3.17) 



3.1.2. Homogeneous SM equation and diffuse transmission 

The solution to the homogeneous SM equation is a priori of the form 



/A h (t) + B h (t)(1-v 2 ) 
Ah(t) 


V o 



(3.18) 



We deduce from the integral equations (3.2) for the functions Ajj{t) and Bjj(t), in the 
absence of source terms, the following equations for their Laplace transforms cih(s) and 
b H (s) 



- two (s) - m 2 (s)^ a H (s) + (m 4 (s) - m 2 (s))b H (s) = A H (s), 
(3m 2 (s) - m (s))a H (s) + - 2m (s) + 5m 2 (s) + 3m 4 (s)^ b H (s) = B H (s), 



(3.19) 



with right-hand sides 
dt 



A H (s) 



Bh(s) = J 



2ni(t - s) 
dt 

2ni(t- s) 



(mo(t) + m 2 (t))a H (t) + (m 2 (t) - m 4 (t))b H (t) 

(m (t) - 3m 2 (t))a H (t) + (2m (t) - 5m 2 (t) - 3m 4 (t))b H (t) 

(3.20) 
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On the other hand, the asymptotic behaviour (2.22) implies 

a H (s) = \-^ + 0(l) ( a ->0), (3.21) 
s z s 

while &#(()) is expected to be finite in this limit. 

The determinant of the 2 x 2 linear system (3.19) can be factorised as (16/9)0i(s)02(s). 
This system can be put in diagonal form by looking for linear combinations of the lines of 
eq. (3.19) involving only <pi(s) or 02 (s) acting on the unknowns. We thus get 

<Pi(s)a H (s) = -^( y (3-2s 2 )A H (s)+B H (s)^ (3.22a) 
Ms) ((1 - s 2 )b H (s) - s 2 a H (s)) = - A (Ah(s) + B H {sj). (3.22b) 

We now solve these equations by means of the so-called Wiener-Hopf technique. We 
consider first eq. (3.22a), and we start by investigating the case where <p\{s) is a rational 
function of the form (3.10), with zeros at s = ±zi ja and poles at s = ±pi,b- We observe 
that an(s) is regular for Res < 0, while the right-hand side of eq. (3.22a) is regular for 
Res > 0. Moreover, eq. (3.20) implies that this right-hand side grows at most linearly as 
s — > — oo. Hence the solution of eq. (3.22a), normalised by the condition (3.21), reads 



JV 

ii(i- s /p hb ) 

1 — cs b=1 1 — cs Hi(s) 

i-r ~ 2 

U(l-s/z ha ) 

a=l 



a H (s) = — 2 = — 2 (3.23) 



where c is a constant, yet to be determined. Similarly, the solution of eq. (3.22b) reads 

b H (s) = ^ ((1 - cs)^ - q H 2 (s^j , (3.24) 

where q is another constant. The notation c and q follows ref. [1]. These two constants are 
determined by expressing that the right-hand side of eq. (3.24) remains finite as s — > ±1. 
We finally obtain 

_ 2HK-1) - H 2 (-l) 2y^g 1 (-l)g 2 (-l) 
^(-l) + 2fl|(-l)' q H*(-l) + 2Hl(-iy 

The representation (3.13), (3.14) permits a numerical evaluation of these numbers, and of 
all the subsequent quantities, with arbitrary accuracy. We thus get 

c = 0.872 941, q = 0.487 827. (3.26) 
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It is worth noticing that the exact solution derived above does not require to determine 
the auxiliary functions Ah(s) and Bh(s) explicitly. 

The observables of interest can now be deduced as follows. 

• The constant tq is obtained by comparing the result (3.23) with the expansion (3.21), 
namely 

To = c - = 0.712110. (3.27) 

v2 

This number is remarkably close to the celebrated value for isotropic scattering of scalar 
waves, recalled in Table 2. 

• The functions Ti(/t) and r 2 {^) are obtained from their definition (2.24), yielding 

tM = |(a H (-l///) + (1 - n 2 )b H (-l/n)), T 2 (fj) = la H (-l/fj,), (3.28) 
i.e., explicitly, 

tM = ^ 2 # 2 (-l///), T2(jjl) = ^=»(» + c^i-l/fi). (3.29) 

In order to make a comparison with the case of scalar waves, we must take into account 
that the above results describe a single polarisation state, and should be compared with 
half the corresponding quantity for isotropic scattering of scalar waves, determined in refs. 
[19, 20], and denoted there by Ti(/t), and hereafter by r sca i(/i). At nearly grazing incidence, 
both functions ri(/t) and T2(/t) vanish linearly, according to 

3 3 
ti(ai) « -g// = 0.731740//, M^) ~ — t=c// = 0.925 893 /*, (3.30) 
2 2y 2 

while in the scalar case we have r sca i( / u)/2 (y / 3/2) / u = 0.866 025 //. At normal incidence, 
both functions take the common value 

ri(1)=T2(1) = |pizW = 2 . 5 3 876l! ( , 31) 

which is again very close to the corresponding number in the case of scalar waves (see 
Table 2). The full functions Ti(/t) and r 2 (/u) are plotted in Figure 1. They hardly differ 
from each other, and from half the corresponding scalar quantity r sca i(/t)/2. 

In order to underline the main novelty with respect to the scalar case [19, 20], namely 
polarisation effects, we plot in Figure 2 the degree of polarisation P, defined in eq. (2.2), 
which reads in the present case 

= 44^44- (3.32) 
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This quantity has a maximum at grazing incidence, namely 

c-qV2 



P(0) 



0.117127, 



(3.33) 



and vanishes at normal incidence, as it should. 

3.1.3. Inhomogeneous SM equation and diffuse reflection 

The special solution of the full inhomogeneous SM equation can be derived by solving 
the six equations (3.2), along the lines of the previous subsection. 

Let us begin with the functions A(r) and B(r). Their Laplace transforms a(s) and 
b(s) still obey equations of the form (3.19), albeit with the contributions of the source 
terms in their right-hand sides: 



dt 



+ 



2jvi(t-s) 



(m (t) + m 2 (t))a(t) + (m 2 (t) - m 4 (t))b{t) 



B(s) = ((2 -3^)I 1 -I 2 ) T 



fJ>a 



(3.34) 



+ 



dt 



2ni{t - s) 



(m (t) - 3m 2 (t))a(t) + (2m (t) - 5m 2 (t) - 3m 4 (t))b(t) 



These equations can be solved by means of the Wiener-Hopf technique, along the lines 
of the previous subsection. The undetermined constants can be fixed in terms of c and q, 
given by eq. (3.25), and we finally obtain 



a(s) = 
b(s) = 



3// ifi(s) 
2y/2s 



rr ( , / \t | + C~ (1 +^ a c)gii"l(-l/^ a ) 

qfjL a H2(-l/fJL a )Li H ; ^ h 



V2 



s 2 a(s) 
S/j, a H 2 (s) 



\l a - C - (1 - H a c)s H^-l/fla) 
— H2(-l/tla)h +Q ^ h 



1 - SH a 



V2 



(3.35) 

The other four functions C(r), • • • , F(t) are easier to determine, since the last four 
lines of eq. (3.2) are uncoupled. We thus obtain the following closed-form expressions for 
their Laplace transforms 



c(s) 
d(s) 
e(s) 

m 



3 fi% 



2 1 - S/l a 

3 fiaU a (2fi a I 1 +il 3 ) 
4 



1 - s/i a 

3 fl a VaU 



H 5 {s)H 5 (-l/[x a ), 



(3.36) 



H^H^-l/n*), 

H 4 (s)H 4 (-l/fi a ). 



4 1 - SfJL a 

3 Haifilh - h + ifiah) 



1 - s/j, a 
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The above results (3.35), (3.36) allow us to give the following expression for the full 
bistatic matrix in the absence of internal reflections: 



7(^a,^a,^6, <Pb) = Yl 7 (fc) (>a,^)e 



ik((p a —<pb) 



with 



k = -2 



(3.37) 



7 (0) (^a,^b) = ^MaM&X 



X 



Mo + Mb 



V2 



l 2 ±± 2 



H$H\ 






Wb H a H b 



V2 

1 + l^a^b 
Ma + Mb 



+ cj HfH\ 

' 






MaMb 
Ma + Mb 



b 

5 -"5 



7 (1) (Ma,Mb) = 7 l lj '(Ma,Mb) = 7 



3 HaVaUbVb 



/-2fi af i b H$Hl 



4 Ha + Hb 





-2in a H%H* 












£7" a o"b 









iff iff / 



7 (2) (^ a ,^ 6 ) =7^ ' 2} ' {jla^b 



(-2)" 



8 /J ia + fib 



( Anl "Mb ^MaMb °\ 

-Ma 1 -«Ma 
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0/ 



(3.38) 



and with the notation H% = H n {— l/fJ, a ), H h n = H n (— 1/m&)- 

At normal incidence (fi a = \i b = 1, <p a = <pb = 0), the above expression simplifies, and 
it agrees with the general form (2.42), with 

ln= T Tt 1) 3Z 1 l + 3 ^ = 3.025 270, 



fl?(-l) + 21f 2 2 (-l) 
3ff 2 (-l)ff 2 2 (-l) 3g|(-l) 

712 iy?(-i) + 2#|(-i) 8 

3if 2 (-l) 
744 = ^7 = -1.086 530. 



1.563100, 



(3.39) 



We now derive a few special results of interest from the above general expressions. 
First, neglecting polarisation effects, the total diffuse reflected intensity in the normal 
direction is given by 
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(3.40) 



This number is slightly above the corresponding value 7(1, 1) for isotropic scattering of 
scalar waves (see Table 2). 

The interesting polarisation dependence of the enhancement factor at the top of the 
backscattering cone, described in general terms in section 2.6, can also be made fully 
quantitative in the present case. For linearly polarised beams, the extremal enhancement 
factors (2.56), corresponding to parallel and perpendicular detection, read 

B\\ = 1.752 088, B± = 1.120 158, (3.41) 

while the width of the triangular cone for parallel detection (2.61) is 

AQ|| = 0.704 063. (3.42) 

For circularly polarised beams, the enhancement factors (2.65) in the helicity-preserving 
channel and in the channel of opposite helicity read, respectively, 

Si = 2, B_ x = 1.250 989, (3.43) 

while the width of the triangular cone in the helicity-preserving channel (2.66) is 

AQi = 0.407487. (3.44) 

The most significant of these numbers are again compared with their analogues for isotropic 
scattering of scalar waves in Table 2. 

3.1.4. Extinction lengths of non-diffusive modes 

The exact solution of the inhomogeneous SM equation in the absence of internal reflec- 
tions, derived in section 3.1.3, also allows us to predict the extinction lengths characterising 
the exponential fall-off of the various non-diffusive polarised components of the intensity 
of radiation, deep in the bulk of a semi-infinite sample. These quantities do not depend at 
all on the index mismatch, so that the results derived below are quite general. 

Let us take for definiteness the example of the component of the intensity described 
by the function D(r), defined in eq. (3.1). Its Laplace transform d(s) is by construction 
analytic for Res < 0. The explicit expression (3.36) shows, however, that it is actually 
analytic in a larger domain, defined by Res < 1//U a and Res < S3, where S3 is the first 
pole of Hs(s), namely the first zero of the kernel function 03 (s). Here first means having 
the smallest real part. We have S3 = 0.914 815 < 1 < 1//U a - The first singularity of d(s) 
is therefore a simple pole at s = S3. We have thus demonstrated the exponential fall-off 
-D(r) ~ exp(— t/£s), with a dimensionless reduced extinction length £3 = I/S3 = 1.093 116. 

More generally, all the extinction lengths are given by the locations of the first singu- 
larities of the corresponding Laplace transforms. We thus obtain 

B(t) ~ E(t) ~ e~ T/£ \ C(r)~e- T/ ^ 5 , D(r) ~ e~ r/l \ F{r) ~ e~ T/e \ (3.45) 
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while the function A(r), pertaining to the diffusive sector, does not fall off, but it rather 
admits the limit value A(+oo) = Ti(/x a )Ii + T 2 (/x a )l2 deep inside a semi-infinite sample. 

Thus there are altogether four different dimensionless extinction lengths, £ n = l/s n 
(n = 1,3,4,5), where s n is the zero of the kernel function 4> n (s) having the smallest 
positive real part. Table 3 gives the exact values of the extinction lengths £ n , together 
with the corresponding approximate values ^ lff , predicted by the diffusion approximation 
[10]. This approximate scheme consists in brutally truncating the kernel functions <f> n (s) 
to the second order of their series expansion in s, given in eq. (3.15). 

3.2. Enhanced backscattering peak 

We now derive analytical expressions for the five functions describing the polarisation 
dependence of the enhanced backscattering cone at normal incidence, according to the 
general formalism exposed in section 2.6, in the absence of internal reflections. We start 
by introducing the following parametrisation 

(A{t) + B{t){1-h*) 
A(t) 


V C(r)fi 

/(iD(r)-E(r))fx\ 




r(°)(r,/*) 



rW(r,^) = -r(- 1 ) (r, /*) = !/ 



V 



D(t) +iE(r) 
iF(r) 



(3.46) 



/ 



T^(r^) = T^\t,/jl) = (G(r)+iH(r)) 



( * \ 

-1 

—2ifx 

V o / 



that slightly differs from eq. (3.1). 

Eq. (2.46) implies that the eight real functions A(t) 7 ■ ■ ■ , H(r) obey the following sets 
of coupled linear integral equations, as shown by the braces 

^ ((Ii + h)e~ T + (M 00 + M 02 ) *A+ (M 02 - M 04 ) * B 

- 2M 13 * D + 2(M 20 + M 22 ) ★ G) , 

- (-(Ii + I 2 )e- T + (M 00 - 3M 02 ) *A+ (2M 00 - 5M 02 + 3M 04 ) ★ B 
+ 2(-2M n + 3M 13 ) *D — 2(M 20 + 3M 22 ) ★ g) , 
(2Mu *A + 2(Mn - Mi 3 ) ★ B 

+ (M 00 + M 02 - 2M 04 + M 20 - 2M 22 ) *D + 2(M n + M 13 + M 31 ) ★ (?) , 
| ((Ii - I 2 )e- T + M 20 * A — M 22 * B 

- (Mn + M 13 + M 31 ) ★ D + (M 00 + 2M 02 + M 04 + M 40 ) ★ g) , 



A = 



B = 



D 



G = 



3 
4 

3 

8 



(3.47a) 
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(l 4 e- T + M 02 ★ C - 2M n * F) , 

V 7 (3.47b) 

(Mh * C + (M 00 - M 02 - M 20 ) * f) , 

((M 00 + M 02 - 2M 04 - M 20 + 2M 22 ) * E + 2(M n + M 13 - M 31 ) * if) , 

(3.47c) 

H = - (l 3 e- T + (-M n - M 13 + M 31 ) * E + (M 00 + 2M 02 + M 04 - M 40 ) * Hj . 
The kernels involved in these equations read 

M pq (Q,r) = I ^^^e-^^J q (Qu\r\/^), (3.48) 
Jo Z H 

where v has been defined in eq. (2.3). 

3.2.1. Preliminaries: kernels and H-functions 

The explicit solution of eq. (3.47) again involves the Laplace transforms a(s), • • • , h(s) 
of the functions introduced in the parametrisation (3.46). In a first step, we must there- 
fore evaluate the Laplace transforms m pq (Q, s) of the kernels M pq (Q, r). Let us take the 
example of moo(Q,s)- The representation (3.48) yields 



m 00 (Q,s)= I ^e(r//)exp ( sr - - + iQ— simp ) , (3.49) 

-i z Jo Zn J-oo 1AM \ A 1 A 1 



ran 

J 4^ 



where © is Heaviside's function, and dfi is the solid-angle element on the unit sphere of 
n, with co-ordinates X = vcofnp = sin#cos</?, Y = z/sin^ = sin6>sin</?, Z = fx = cos 9. 
Performing the r-integral yields 

m 00 (Q,s) = [ ^- 1 . (3.50) 

J 4% 1 — s cos v — iQ sin 6 cos <p 

The denominator can be transformed by a suitable rotation in the Y-Z plane into 1 — aZ' , 
with 

a = ^s 2 -Q 2 . (3.51) 

We thus obtain 

moo(Q, s) = m (a). (3.52) 
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It turns out that all the kernels involved in eq. (3.47) can be treated similarly, and expressed 
as linear combinations of the m 2p (cr), which have been determined in eq. (3.7), namely 



m 00 (Q, s) = m (a). 



Q 2 



mu(Q, s) = ^(3m 2 (a) - m (a)). 



m 02 {Q, s) = m 2 (cr) + ^(3m 2 (cr) - m (cr)), 

Q 2 Q 4 
m 4(Q, s) = m 4 (cr) + ^-(5m 4 ((r) - 3m 2 (a)) + ^-j (35m 4 (cr) - 30m 2 (a) + 3m (a)). 

Qs 

2(7 S 

Qs Q 3 s 
mis(Q, s) = 7^2 (5m 4 (cr) - 3m 2 (cr)) + -^j (35m 4 (cr) - 30m 2 (a) + 3m (cr)), 

Q 2 

m 20 (Q, s) = ^(3m 2 (cr) - m (a)), 

Q 2 <5 4 

m 2 2(Q, s) = ^-^(l5m 4 (cr) - 12m 2 (a) + m (a)) + ^j(35m 4 (a) - 30m 2 (cr) + 3m (a)), 
Q 3 s 

"7-31 {Q, s) = -^j(35m 4 (a) - 30m 2 (a) + 3ra (cr)), 
^4o(<3, s) = — ^ (35m 4 (<7) - 30m 2 (a) + 3m (a)). 

(3.53) 

We are again led to consider the kernel functions <j> n (cr), defined in eq (3.8). Although 
the latter only depend on the variable a of eq. (3.51), the associated .//-functions H n (Q, s) 
depend separately on both variables Q and s. The latter functions are still defined by the 
Wiener-Hopf identity (3.9), which now reads 

*" (ff) = Hn (Q.s)H n (Q.-sy (3 ' 54) 

with the condition that the H n (Q,s) be analytic in the left half-plane Res < 0. For a 
rational kernel function <f>(cr) of the form (3.10), the .ff -function reads 

N 



6=1 V ' 



H(Q, s) = ^ . (3.55) 



M 



a=l 



In the present case, the functions H n {Q, s) still possess the integral representation (3.13), 
up to the replacement 



(3 -> 0(Q, /3) = arctan sj Q 2 + tan 2 13. (3.56) 
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3.2.2. Derivation of 744 (Q) 

We begin with the exact solution of the 2x2 system (3.47b), yielding, after a Laplace 
transformation, c(Q,s) and f(Q,s), and 744(Q) through c(Q, — 1) = — 744(^)14. The 
determinant of this linear system can be factorised as — (9/8)0i(cr)05(cr). Along the lines 
of section 3.1, the system is made diagonal by setting 



c(Q, S ) = sc(Q,s)-2Qf(Q, S ), 
f(Q,s) = (Q/2)c(Q,s)-sf(Q,s). 

The inverse formulae read 

a 2 c(Q,s) = sc(Q,s) - 2Qf(Q,s), 
a 2 f{Q,s) = {Q/2)Z{Q,s)-sf{Q,s), 

and the new functions obey 

</> 5 (a)c(Q, s) = z ~ I4 + C(Q, s), 
1 — s 

1 — s 

where C{Q, s) and ^F(Q, s) are defined in analogy with eq. (3.20). 
The solution of eq. (3.59) reads 



(3.57) 



(3.58) 



(3.59) 



c(Q,s)= (ci + ^^j UH 5 {Q,s), 
f(Q,s)=(fi + ^A UH x {Q,s). 



(3.60) 



This expression involves, besides the corresponding if-functions, four constants, yet to be 
determined. The s — > 1 limit of eq. (3.59) fixes two of them: 

c 2 = ^UH 5 (Q, -1), h = iQUH^Q, -1). (3.61) 

The last two constants, c\ and /1, are then fixed by requiring that c(Q,s) and f(Q,s), 
as given by eq. (3.58), are finite for a — > 0, i.e., s Q and s — > —Q. Skipping lengthy 
details, we finally get 

744(g) " 4(l-Q^(HHQ,-Q) + Hl(Q,-Q))> (3 - 62) 

where 

AT 44 (Q) = Q 2 (l + Q) 2 Hl(Q, -Q)H 2 (Q, -1) + Q 2 (l - Q) 2 H 2 (Q, -Q)H 2 (Q, -1) 

+ (1 - Q) 2 H 2 (Q, -Q)H 2 (Q, -1) + (1 + Q) 2 H 2 (Q, -Q)H 2 (Q, -1) (3.63) 
- 8Q 2 H 1 (Q, -Q)H 5 (Q, -Q)#i(Q, -1)# 5 (<2, -1). 
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3.2.3. Derivation of 733 (Q) 

The 2x2 linear system (3.47c) yields, after a Laplace transformation, e(Q, s) and 
h(Q,s), and 733 (Q) through h(Q, — 1) = (1/4)733 (Q)lz- The determinant of this system 
can be factorised as (9/32)03 (cr)04(cr). Its exact solution closely follows the lines of the 
previous subsection. We are led to consider the linear combinations 

e(Q, s) = -se(Q, s) + 2Qh(Q, s), 
HQ, s) = -(Q/2)e(Q, s) + sh(Q, s), 

which obey 

03(<T)e(Q,s) = --^-I 3 + 5(Q,s), 

(3.65) 

M<?)HQ, s) = ~r~ 1 3 + H(Q, s). 
1 — s 

The solution of these equations is fully similar to that of eq. (3.59). We are thus left with 

r n \_ _3 A/" 33 (Q) . . 

733WJ " *(i-Q 2 ) 2 (Hi(Q,-Q) + Hi(Q,-Q)y [6 - bb) 

and with 

AT 33 (Q) = Q 2 (l + Q) 2 Hj(Q, -Q)H 2 3 (Q, -1) +Q 2 (1- Q) 2 Hl(Q, -Q)H 2 (Q, -1) 

+ (1 - Q) 2 H 2 (Q, -Q)H 2 (Q, -1) + (1 + Q) 2 H 2 (Q, -Q)H 2 (Q, -1) (3.67) 
- 8Q 2 H 3 (Q, -Q)H 4 (Q, -Q)H 3 (Q, -1)H 4 (Q, -1). 



3.2.4. Derivation of 711 (Q), 712 (Q), and 722 (<5) 

The determinant of the 4x4 linear system (3.47a), after a Laplace transformation, 
can be factorised as (81/256)0i(u)02(c)03(c)04(c)- Its exact solution, along the lines of 
the previous cases, involves algebraic manipulations on very lengthy expressions. Some of 
them have been either carried out, or just checked, by means of the MACSYMA software. 
Furthermore, the final step of this calculation, involving the solution of the 9x9 linear 
system (3.74), must for practical purposes be performed numerically. 

In a first step, the system (3.47a) is put in diagonal form by introducing the linear 
combinations 

a(Q, s) = 2a 2 a(Q, s) - Q 2 b(Q, s) - 2Qsd(Q, s) + 2Q 2 g(Q, s), 
b(Q, s) = -2a 4 a(Q, s) + (-2a 4 - 2Q 2 a 2 + 2a 2 + 3Q 2 )b(Q, s) 

+ 2Qs(-2a 2 + 3)d(Q, s) + 2Q 2 (2a 2 - 3)g(Q, s), (3.68) 
d(Q, s) = Qsb(Q, s) + (a 2 + 2Q 2 )d(Q, s) - 2Qsg(Q, s), 
9(Q, s) = (Q 2 /2)b(Q, s) + Qsd(Q, s) - (2a 2 + Q 2 )g(Q, s). 
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These new functions obey 



rh < \~(n \ 2( g 2 -l)(I 1+ I 2 )-2Q 2 I 2 ~ 

0i(t7)o(<5, s) = + A(Q, s), 

1 — s 



^(<T)b(Q,s)='^^ + B(Q,s), 



1 T ( 3 - 69 ) 

<j> 3 (a)d(Q t s) = -^-j + V(Q,s), 

\~(n \ 2 g 2 (I 1 -I 2 )+2Q 2 I 2 ~ 

J. S 

where ^4.(Q, s), £>(<2, s), ^(<5, s), and Q{Q, s) are defined in analogy with eq. (3.20). 
The formulae inverse to eq. (3.68) read 

4a 4 ((T 2 - l)a(Q, s) = a 2 (2a 2 + Q 2 - 2)a(Q, s) + Q 2 b(Q, s) 

+ 4Qs(a 2 - l)d(Q, s) - 2Q 2 (a 2 - l)g(Q, s), 
4a 4 ((T 2 - l)b(Q, s) = -a 2 (2a 2 + 3Q 2 )a(Q, s) - (2a 2 + 3Q 2 )b(Q, s) 

- 12Qs(a 2 - l)d(Q, s) + 6Q 2 (a 2 - l)g(Q, s), 
2a 4 (a 2 - l)d(Q, s) = Qsa 2 a(Q, s) + Qsb(Q, s) 

+ 2(a 2 - l)(a 2 + 2Q 2 )d(Q, s) - 2Qs(a 2 - l)g(Q, s), 
8a\a 2 - l)g{Q, s) = Q 2 a 2 a(Q, s) + Q 2 b(Q, s) 

+ AQs{a 2 - l)d(Q, s) - 2{a 2 - l){2a 2 + Q 2 )g(Q, s), 

so that we have 



7n (0)1, + 7o <«)!, = ~ b -^ + iQd(Q ' + 2diQ ' 



72l(Q)Il +722(<5)I 2 = 



2(1 - Q 2 ) 2 
a(Q,-l) + 2g(Q,-l) 



2(1 - Q 2 ) 
The solution to eq. (3.69) reads 



a(Q. .s) = ( ai + a 2 s + Hx(Q, s), 

b(Q, s)=[b 1 + b 2 s + b 3 s 2 + ^-^j H 2 (Q, s), 
d{Q, s)=[d 1 + d 2 s + H 3 (Q, s), 

g(Q, s) = [ gi + g 2 s + j H 4 (Q, s), 
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(3.70) 



(3.71) 



(3.72) 



where a±, ■ ■ ■ , gs are 13 Q-dependent constants to be determined. 
The s — > 1 limit of eq. (3.69) fixes four of these constants: 

a 3 = -fahH^Q, -1), b 4 = -3Q 2 hH 2 (Q, -1), 

3 3 (3-73) 

rf 3 = --QIiH s (Q, -1), ^3 = - i (li + (Q 2 - l)h)H 4 (Q, -1). 

The last nine constants are then determined by expressing that the functions a(<5, s), • • •, 
g{Q, s) have the expected regularity properties at the points where the inversion formulas 
(3.70) are singular, namely a 2 = 0, i.e., s = ±Q, or a 2 = 1, i.e., s = ±a/1 + Q 2 . We thus 
obtain the following system of nine linear equations 



a (Q, -y/l + Q 2 J+b (Q, - Vl + Q 2 ) = 0, 


(3.74a) 


a (Q, v 7 ! + Q 2 )+b (Q, v 7 ! + Q 2 ) = 0, 


(3.74b) 


a(Q,-Q) + 2g(Q,-Q)=0, 


(3.74c) 


KQ, -Q) - WQ, -Q) = o, 


(3.74d) 


-Q) + 2^(Q, -Q) = 0, 


(3.74e) 


(6(Q, s) + 4d(Q, s) + 2g(Q, s)) - 8Qg(Q, -Q) = 0, 
as V / S =-Q 


(3.74f) 


6(Q,Q)-6^(Q,Q) = 0, 


(3.74g) 


d(Q,Q)-2g(Q,Q) = 0, 


(3.74h) 



3^-f6(Q,s)-4d(g,s) + 2^(Q,s)) + 10Qa(Q, Q) + UQg(Q, Q) = 0, (3.74i) 

as V / s=q 

where (a) and (b) express the regularity of the functions a(Q, s), ■ ■ ■ , g(Q, s) (absence of 
pole) at s = ±a/1 + Q 2 ; (c) to (f) express their regularity (absence of double and of simple 
pole) at s = —Q; (g) to (i) express the absence of double pole at s = Q, as well as the 
proportionality of the residues of the simple poles to the right null vector of the system 
(3.47a), V = (3 - 2Q 2 , 3 + Q 2 , 2Q 2 , Q 2 /2). 

The last of the above properties implies that the four functions A(r), B(t), D(t), and 
G(t) fall off as exp(— Qr). The dimensionless extinction length in the diffusive sector thus 
reads £(Q) = 1/Q in units of the mean free path £, i.e., 

m = ^ = \ (3.75) 

in physical units. This simple result holds for any value of the transverse wavevector q. 
Moreover, since it is a bulk property of the problem, it also holds in the presence of an index 
mismatch, just as the extinction lengths of the non-diffusive sectors for Q = 0, determined 
in section 3.1.4. 
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By inserting the explicit forms (3.72) into the expressions (3.74), we obtain a 9 x 9 
linear system for the £?-dependent constants {ai, a2, 6i, 62, 63, g^l, ^2, <7i, <72}, whose coef- 
ficients have complicated expressions involving the functions H(Q,s). This system has 
been solved formally by means of the MACSYMA software: the outcome for each constant 
contains thousands of products of up to seven ii/"-functions, so that this approach is of no 
practical use. The above system is however easily solved numerically, for any given value 
of Q. This is the way we have chosen to follow for practical purposes. 

3.2.5. Summary of results 

We have achieved the exact analytical determination of the enhanced backscattering 
cone in the absence of internal reflections. Its dependence on polarisations is contained in 
five functions of the reduced wavevector Q. Two of them, 733(6?) and 744(6?), are given 
explicitly in eqs. (3.62), (3.66), while the other three, 711(6?), 712 (6?), and 722 (6?), are 
determined analytically, up to a last step which consists in solving numerically a well- 
posed 9x9 linear system, for any fixed value of 6?. 

The following regimes are of special interest. 

• For small Q, i.e., in the vicinity of the top of the cone, 711(6?), 712(6?)) and 722(6?) have 
the common characteristic triangular shape given by the general result (2.58), while 733(6?) 
and 744(6?) have a smooth dependence on Q 2 . The enhancement factors and widths of 
the triangular cone, for linear and circular polarisations, have already been given in eqs. 
(3.41-44). 

• For large Q, i.e., in the wings of the cone, the leading contributions come from low-order 
scattering events, as usual. The single-scattering contribution (2.45) is of little interest, 
since it is subtracted in the formula (2.41) for the enhanced backscattering peak. The 
leading large-C? behaviour of the enhancement factor is thus given by double-scattering 
events. The contribution of this class of events can be obtained by solving eq. (3.47) to 
first order in the kernels M pq . As it turns out, for large Q these kernels become small, as 
expected, but also local in the r-variable: M pq (r, r') ~ m pq (Q,0)5(T — r'). Furthermore, 
only the following kernels contribute to leading order in 1/6?: 

7T 7T 37T 7T 7T 3lT 



moo "2Q' mo2 "4Q' mo4 "l6Q' m20 "4Q' m22 "l6Q' ^ ~ 166?' 

(3.76) 

so that we are left with 

(3.77) 

We end up by illustrating a few interesting features of our results. We first consider 
linearly polarised beams at normal incidence and for parallel detection, in the following 
two geometries: 

(i) both polarisations parallel to Q, i.e., ip a = ipt, = 0, 

(ii) both polarisations perpendicular to Q, i.e., ip a = iftb = tt/2. 
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The enhancement factors, given by eq. (2.53), namely 



B « (Q) = 1 + lliM__^, BjM ) (0) = 1 + ^(0)-3/4 
ii 7m u ) 11 7iH u ) 

are plotted in Figure 3. Both curves coincide with the result (3.41) at Q = 0, while they are 
slightly different from each other at Q ^ 0. The enhancement factor of isotropic scattering 
of scalar waves, after refs. [30, 19], is also shown for comparison. 

Similarly, we consider linearly polarised beams with perpendicular detection, in the 
following two cases: 

(iii) one polarisation parallel to Q, i.e., ip a = 0, ipb = n/2, 

(iv) both polarisations at 45° with respect to Q, i.e., ip a = —ipb = tt/4. 
The enhancement factors, 

u( m hn\ — i _l 744(Q) — 733(Q) c>( iv )/m _ -\ , T 11 

(Q) + 72 2 (Q) - 2712(g) + 2 744 (Q) 

^ {Q) ~ 1+ 2^(0) ' ^ 4^ > 

(3.79) 

are plotted in Figure 4. Both factors coincide with the result (3.41) at Q = 0, while they 
are slightly different from each other at Q 7^ 0. 

We end up by considering circularly polarised beams at normal incidence. The en- 
hancement factors Bi(Q) of the helicity-preserving channel, and B-i(Q) of the channel of 
opposite helicity, given by inserting the above results into the general expressions (2.54), 
(2.63), are plotted in Figures 5 and 6, respectively. 

4. LARGE INDEX MISMATCH REGIME 

Previous works [19-21] on multiple scattering of scalar waves suggest that the SM 
equations cannot be solved analytically in the presence of internal reflections, which take 
place whenever the index mismatch m = njn\ is different from unity. However, and 
interestingly enough, it has been shown in refs. [19-21] that the problem becomes again 
tractable by analytical means in the regime of a large index mismatch (m<lorm>l), 
at least in the case of scalar waves. The intuitive origin of this simplification is as follows. 
If the ratio m of both optical indices is very small (respectively, very large), there is 
total reflection for almost any incidence angle outside the medium (respectively, inside the 
medium), except for a narrow cone around the normal incidence. As a consequence, the 
radiation undergoes many internal reflections at the boundary of the sample, and hence 
many scattering events, before it can exit the medium. In this section we extend this 
approach to the Rayleigh scattering of electromagnetic waves. 

4.1. Diffuse reflection and transmission 

In this section we investigate the diffuse intensity in reflection and in transmission 
in the large index mismatch regime. To do so, it is advantageous to consider the matrix 
Green's function G(r, (p, r', //, </?'), which obeys eq. (2.13). We anticipate on physical 
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grounds that only the ^-independent sector (k = 0) is important. We thus rewrite eq. 
(2.13) as 

Jo Jo Z A* 

+ ir ^ Jo |^ e " (T "" T)/ ' i " P(0)( ^ VO-G^Cr", V, r', „') 

+ / + °°dr" / 1 ^ e -^")/,"p(0) (ft/) . G (0) (/)/)T>1 
JO Jo ^ 

- f +0 ° ^" ^e-^")/-" (1 - R(^'))-P (0) (^ /).G(°)(r", r', „'). 
Jo Jo 

(4.1) 

The above observations suggest to treat the small matrix 1 — R(/u) as a perturbation. 
In the limit of an infinite index mismatch (m = or m = +oo), this matrix vanishes 
identically. The rest of eq. (4.1), without the last line, has a zero mode of the form 



M nP , = I 



nat A nat ^ -"-nat 



Irmt 



/l 1 0\ 

110 



Vo o o 0/ 



(4.2) 



independent of r and [i. This is demonstrated by the identity 

11 d/i. 



J ^ ^M nat .p(°)(/i, n') = M nat for all n' . 



(4.3) 



Along the lines of refs. [19-21], we expect that the Green's function becomes propor- 
tional to the constant matrix M nat , with a diverging prefactor, in the large index mismatch 
regime. We thus look for a singular expansion of the form 



G<°>(r, (jl, t', fj,') = CM nat + G< 0) (r, ^ r', »') + 



(4.4) 



(0) 



where it is understood that the constant C diverges as m — * or m — > +oo, while Gq 
stays finite, and the dots stand for higher-order corrections. We insert this expansion into 
eq. (4.1), and then act on both sides with the operator j + °° dr J^ 1 (d / u/2)M nat . Integrals 

over the finite part of the Green's function cancel out, so that we are left with a 

simple expression for the constant C, namely 



with 



T 



77, 



f 

Jo 



2//d//Tii(//), 71 



(4.5) 



(4.6) 
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The above quantities only depend on the index mismatch m. They are interpreted as the 
mean flux transmission coefficients of one boundary of the sample, averaged over incidence 
angles. 7j| and 71 correspond to prescribed polarisations, while 7" is also averaged over 
both polarisation channels. 

More explicitly, the expressions (2.8), (2.9) of the Fresnel intensity coefficients allow 
us to perform the integrals (4.6) in closed form, for both m > 1 and m < 1. It turns out 
that both cases can be gathered in the following formulas, valid for m > 1: 

4(2m+ 1) 
3m(m + l) 2 ' 
2m(m 2 - l) 2 , m + 1 

— t—z r~ — In (4.7) 

(m 2 + l) 3 m-1 y 1 

16m 3 (m 4 + l) , 4m 2 (m 2 + 2m — 1) 

lnm + 







m 


\m 






m 1 


V m 



(m 2 -l) 2 (m 2 + l) 3 (m 2 -l)(m 2 + l) 2 ' 

As expected, the flux transmission coefficients vanish in the regime of a large index mis- 
match, according to 

. _. 8m _ 16m 

m -> : 7j| sa 8m, 71 « — , T « 

^ 8 8 ^ 16 ( 4 - 8 ) 

m — > +oo : Tii ps — ? , Tj_ T 

m 



3 ' ^'^ 3m 3 ' 3m 3 ' 



The knowledge of the matrix Green's function yields by eqs. (2.16), (2.22-24) the 
following predictions for observables of interest in the large index mismatch regime 

^W)^ (U = l,2), 

riO*)«^ (< = 1, 2), (4.9) 
4 

T0 "3T- 

These leading-order results in the large index mismatch regime are very similar to those 
obtained in the case of isotropic [19, 20] and arbitrary anisotropic [21] scattering of scalar 
waves. Figure 7 shows plots of the mean transmission amplitudes T, 7j|, and Tj_, against 
the index mismatch m. It is worth noticing that the reflection and transmission coefficients 
for scalar waves coincide with R±((jl) and T±(fj,), so that 71 was already involved in the 
predictions of refs. [19-21] for scalar waves. 

The behaviour of the quantities investigated above in the large index mismatch regime 
involves, besides the leading asymptotic behaviour in 1/7" derived above, finite parts re- 
lated to the Green's function G'q K Moreover, also governs the (non- divergent) large 
index mismatch behaviour of all the other quantities, like e.g. the entries of the bistatic 
matrix 7^ outside the (1,2) sector, or the bistatic matrices 7^ for k 7^ 0. These finite 
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parts cannot be determined analytically in general. In the case of isotropic scattering of 
scalar waves, the finite parts of some simple observables have been determined [19, 20], 
either analytically or numerically. 

4.2. Enhanced backscattering cone 

We now extend the above analysis to the enhanced backscattering cone. In analogy 
with the case of scalar waves, treated in refs. [19-21], we want to show that the bistatic 
matrix 'j(Q) takes a simple scaling form in the regime of small Q and large index mismatch, 
where enhanced backscattering is dominated by long-distance effects. 

To do so, we look for a solution to the Q-dependent SM equations (2.46) in the form 

T« (r, » a(Q)e- Qr M nat 4o, (4.10) 

for m < 1 or m > 1, and Q <C 1. The above ansatz is justified as follows: the exponential 
fall-off in exp(— Qr) is quite general [see eq. (3.75)]; the ^-dependent sectors (k ^ 0) 
are again neglected; the proportionality to the matrix M nat is assumed because of the 
structure (4.2), (4.3) of the zero mode of the RTT problem at Q = 0. 

In analogy with section 4.1, we insert the form (4.10) into eq. (2.46), and then act on 
both sides with the operator J + °° dr J^ 1 (d / u/2)M nat . The integrals which do not involve 
the small matrix 1 — R(/u) can be performed exactly; their Q-dependence is given by 
woo(<3,0) = mo(iQ) = arctan(Q)/(5 ~ 1 — Q 2 /3, by eq. (3.7). The Q-dependence of the 
integral involving 1 — R(/i) can be neglected. Consistently neglecting all corrections of 
order Q 2 , we obtain 

liAQ) ~ 2Q l T (<,J = 1,2). (4.11) 
~Y + ~2 

This prediction of the leading behaviour of the enhanced backscattering cone in the large 
index mismatch regime is again very similar to the case of scalar waves [19-21]. 

We end up by giving a few consequences of the above predictions in the large index 
mismatch regime. In the case of linear polarisations, the enhancement factor (2.55) at the 
top of the cone reads B(0) = 2 cos 2 It assumes the maximal value B\\ = 2 for parallel 
detection, because the single-scattering contribution is negligible. The width of the top of 
the cone (2.61) reads 




(4.12) 



This last result also gives the width AQi of the top of the cone of enhanced backscattering 
for the helicity-preserving channel (2.66) in the case of circular polarisations. 
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5. DISCUSSION 



In this paper we have extended to the Rayleigh scattering of electromagnetic waves 
our previous investigations [19-21] of multiple scattering of scalar waves in a thick-slab 
geometry. Both the setup and the formalism of the present work closely follow those refer- 
ences, so that only the main lines of the derivations have been reproduced here. The main 
advantage of this approach, based on RTT, is that the role of skin layers, and especially 
the effects of internal reflections, are incorporated in a natural way. The present approach 
has no phenomenological or approximate character, besides the restriction of its validity 
to the regime A <C I <C L, in contrast with the widely used diffusion approximation. Only 
few analytical results [17, 18] had been obtained for electromagnetic waves along this line 
of thought, since the pioneering work of Chandrasekhar [1]. It is, however, worth mention- 
ing that the vector RTT formalism, including the effects of internal reflections, has been 
exposed earlier [31, 32], although these authors only solved the SM equations numerically 
in some specific situations, rather than investigating their general properties. We have first 
derived general results on vector RTT in section 2, where the mean values of observables 
are expressed in terms of solutions to vector SM equations, including the effects of polar- 
isations and of internal reflections. Closed-form expressions for these general predictions 
are then derived in two cases, namely in the absence of internal reflections (in section 3), 
and in the regime of a large index mismatch (in section 4). 

In the absence of internal reflections (m = n/ni = 1), the SM equations have been 
solved by means of the Wiener-Hopf technique. We have presented in section 3 a self- 
contained exposition of all the exact results known so far [1] . More importantly, we have 
given the first complete analytical derivation of the cone of enhanced backscattering, com- 
pleting thus the analytical results of refs. [17, 18], as well as some less accurate estimates, 
obtained either by means of the diffusion approximation or by numerical simulations [10, 
11, 13]. As a general rule, illustrated by the first three items of Table 2, quantities which 
are either averaged over the polarisation degrees of freedom, or do not depend on polarisa- 
tions at all, are found to be very close to the corresponding figures in the case of multiple 
isotropic scattering of scalar waves. A similar observation has been made in ref. [21], where 
various observables were compared for isotropic and very anisotropic (forward) scattering 
of scalar waves. The prototype of such quantities is the thickness tq of a skin layer, ex- 
pressed in units of the transport mean free path £*. This number is hardly sensitive to the 
anisotropy of the scattering mechanism nor to polarisations: it always comes out to read 
r « 0.71 [4]. 

Other observables, such as the detailed shape of the cone of enhanced backscattering, 
have a more or less pronounced dependence on the polarisation channels of the incident 
and detected beams. The last two items of Table 2 illustrate this point. The first quantity 
under consideration is the maximal enhancement factor, right at the top of the cone. The 
value S|| of eq. (3.41), corresponding to linear polarisations and parallel detection, as 
well as the value Bi of eq. (3.43), corresponding to circular polarisations and detection 
in the channel of same helicity, are compared to the analogous result for scalar waves 
with isotropic scattering [30, 19], denoted by B: the figures are definitely different from 
each other, although relatives differences are less than 10%. Second, the width of the 
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triangular top of the backscattering cone is considered. The value AQy of eq. (3.42), 
corresponding to linear polarisations and parallel detection, as well as the value AQi of 
eq. (3.44), corresponding to circular polarisations and detection in the channel of same 
helicity, are compared to the analogous result for scalar waves with isotropic scattering 
[30, 19], denoted by AQ: relative differences are more important in this case, going up to 
some 40%. Finally, so far there are essentially no analytical results concerning the RTT 
approach to the general problem of multiple scattering of electromagnetic waves, taking 
into account the combined effects of anisotropic scattering and polarisations. We can infer 
from the results of ref. [21] on the multiple scattering of scalar waves that the anisotropy 
of the scattering mechanism will have little residual effects, once the principal scaling is 
taken into account by expressing observables in terms of the transport mean free path £*. 

In the presence of internal reflections (m = njn\ ^ 1), analytical predictions for 
the various observables of interest have been derived in the large index mismatch regime 
(m C 1 or m > 1), along the lines of previous investigations of scalar waves, with isotropic 
[19, 20] and arbitrary anisotropic [21] scattering. The effects of internal reflections have 
been studied [22-25] using several variants of the diffusion approximation. Ref. [26] 
provides a recent overview of these approaches to the subject. Within the framework 
of RTT, the drastic simplification which takes place in the large index mismatch regime 
has a clear intuitive explanation. Since the transmission through the boundaries of the 
sample is small, radiation is reinjected many times before it can leave the medium. As a 
consequence, the skin layers become very thick and, more importantly, the radiation field 
is uniform over most of these layers. The results (4.9) turn out to have the very same form 
as for multiple scattering of scalar waves, either with isotropic [19, 20] or very anisotropic 
[21] scattering. Most certainly, the very same analytical forms also hold true in the more 
general case of multiple scattering of electromagnetic waves, including both anisotropy and 
polarisation effects, and they are expected to provide an overall satisfactory description of 
the full dependence of physical quantities on the index mismatch, especially in the range 
of most interest (m > 1). 

The present investigations of multiple Rayleigh scattering of electromagnetic waves, 
as well as the previous ones on isotropic and anisotropic scattering of scalar waves [19-21], 
have so far only dealt with the mean intensity, averaged over the random positions of 
the scatterers in the sample. For any given sample of scattering medium, however, the 
intensity has strong fluctuations, which manifest themselves as speckles. For instance, the 
probability law of the fluctuating intensity at a given point is known as Rayleigh's law: 
p(I) ~ exp(— //(/)). The generalisation of Rayleigh's law to polarised radiation is fully 
characterised by the four Stokes parameters [33]. Various correlation functions, aiming 
at a more detailed description of intensity fluctuations and speckle patterns, have been 
the subject of recent theoretical and experimental investigations [34]. We mention the 
extension of these results, in order to include polarisation effects, as an interesting open 
problem. 
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CAPTIONS FOR TABLES AND FIGURES 

Table 1: Definitions and notations for kinematic and other useful quantities. 

Table 2: Comparison of various quantities of interest, denned in section 2, from the 
known exact solutions in the absence of internal reflections. First row: isotropic scattering 
of scalar waves, after ref. [19]. Second row: Rayleigh scattering of electromagnetic waves 
(section 3 of this work). Third row: relative difference of second case with respect to first 
one. 

Table 3: Dimensionless reduced extinction lengths of the various polarised components 
of the diffuse intensity. First row: exact extinction lengths, deduced in section 3.1.4 from 
the solution to the SM equation. Second row: approximate extinction lengths, obtained 
by means of the diffusion approximation [10]. Third row: relative difference of second case 
with respect to first one. Fourth row: notations used in ref. [10]. 

Figure 1: Plot of angular dependence of transmitted intensity in the absence of internal 
reflections, against fx = cos 9. Full lines: ti(/j) (lower curve) and T2(//) (upper curve), 
corresponding to both polarisations channels for Rayleigh scattering (this work). Dashed 
line: r sca i(//)/2, corresponding to half the result for isotropic scattering of scalar waves, 
after ref. [19]. 

Figure 2: Plot of angular dependence of degree of polarisation P(p) of diffuse transmitted 
intensity, in the absence of internal reflections, against fi = cos 9. 

Figure 3: Plot of enhancement factors B\\(Q) for linearly polarised beams at normal 
incidence and parallel detection, in the absence of internal reflections, in two cases defined 
in the text. Upper full line: case (i) (ifj a = ^6 = 0). Lower full line: case (ii) (ip a = ipb = 
7r/2). Dashed line: same quantity for isotropic scattering of scalar waves, after ref. [19]. 

Figure 4: Plot of enhancement factors B±(Q) for linearly polarised beams at normal 
incidence and perpendicular detection, in the absence of internal reflections, in two cases 
defined in the text. Lower full line: case (iii) (ip a = 0, ipb = tt/2). Upper full line: case 

(iv) (V>« = -fa = 7T/4). 

Figure 5: Plot (full line) of enhancement factor B\(Q) for circularly polarised beams at 
normal incidence in the helicity-preserving channel, in the absence of internal reflections. 
Dashed line: same quantity for isotropic scattering of scalar waves, after ref. [19]. 

Figure 6: Plot of enhancement factor B-i(Q) for circularly polarised beams at normal 
incidence in the channel of opposite helicity, in the absence of internal reflections. 

Figure 7: Plot of mean flux transmission coefficients, against optical index mismatch m. 
Upper dashed line: T\\{m). Lower dashed line: T±(m) (also corresponds to scalar waves). 
Full line: their average T(m). 
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Table 1 





outside medium 


inside medium 


optical index 


ni 


n = mri\ 


wavenumber 


ki = niu/c = 27r/Ai 


k = nuj/c = 2/t/A 


incidence angle 


cos 6*i = \J\ — m 2 u 2 

Q1T1 Hi Till/ 

Ollll/l 1 1 LIS 




COS = fX 

c\ri f) 77 a /l //^ 

alii U Is \l J. LI 


parallel wavevector 


p = ki cos 0i 


P = k cos 6* 


total reflection 
condition 


m < 1 and sin6>i > m 
(i.e. P imaginary) 


m > 1 and sin 6> > 1/m 
(i.e. p imaginary) 



transverse wavevector 


q = q = k\ sin 0\ = k sin = kv 


azimuthal angle 
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Table 2 





isotropic scattering 
(scalar waves) 


Rayleigh scattering 
(electromagnetic waves) 


A(%) 


skin layer 
thickness 


r = 0.710446 


r = 0.712110 


0.23 


diffuse 
transmission 


|r Bca i(l) = 2.518 237 


Ti(l) = 2.538 761 (»=i,2) 


0.81 


diffuse 
reflection 


7(1,1) = 4.227681 


7n + 7i2 = 4.588 369 


8.5 


enhancement 
at top of cone 


B = 1.881 732 


Si = 2 


6.3 


width 
of top of cone 


AQ = 1/2 


AQ|| = 0.704 063 
AQi = 0.407487 


41 
-19 
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Table 3 



exact 


diffusion approximation 


A(%) 


notation 
ofref. [10] 


h = l 


if ft = = 0.447 213 


-55 




£3 = 1.093116 


/ 13 

if ft = J — = 0.786 796 
3 V 21 


-28 




£ A = 1.349 587 


/ 93 

if ft = J — = 1.046 536 
4 V 21 


-22 


£2 


£ 5 = 1.172 669 


£f ft = y| = 0.774 596 


-34 


£3 
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